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ABSTRACT 


The first part of this report presents a collection of typical 3 -body . 
trajectories from the libration point on the Sun-Earth line to the Earth. 
These trajectories in the Sun-Earth system may be grouped into four 
distinct families which differ in the transfer time and AV requirements. 
Also included are curves showing the variations of AV with respect to 
transfer time and typical 2 and 3-impulse primer vector histories. 

The second part of the report deals with the development of a 4-body 
trajectory optimization program to compute fuel optimal trajectories 
between the Earth and a point in the Sun-Earth-Moon system. It presents 
methods for generating fuel optimal 2 -impulse trajectories which may 
originate at the Earth or a point in space and fuel optimal 3-impulse 
trajectories between 2 points in space. The extrapolation of the state 
vector and the computation of the state transition matrix are accomplished 
by the Stumpff-Weiss method. The cost and constraint gradients ace 
computed analytically in terms of the terminal state and the state transition 
matrix. The 4-body Lambert problem is solved by using the Newton- 
Raphson method. An accelerated gradient projection method is used to 
optimize a 2-impulse trajectory with terminal constraint. The Davidon's 
Variance Method is used both in the accelerated gradient projection method 
and the outer loop of a 3-impulse trajectory optimization problem. This 
method is preferred over many others mainly because it does not require 
a one -dimensional search. Several well-known methods which have been 
successful in solving 2 -body trajectory optimization problems perform 
poorly in the 4-body system. A brief qualitative comparison of these 
methods is given. 

An example of a 4-body 2 -impulse transfer from the libration 
point to the Earth is included. The difference between this trajectory 
and a 3-body trajectory of the same transfer is readily discernable. 
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INTRODUCTION 


This report is divided into two parts. The first part is a collection 

of typical 3 -body trajectories from the libration point on the Sun-Earth 

line to the earth. They are generated using routines based on the program 

developed by D'Amario^. The extrapolation of the state vector is 

(8 9) 

accomplished by the Wilson's version of the multi-conic method' ’ 

These trajectories in the Sun-Earth system may be grouped into four 
distinct families which differ in the transfer time and AV requirement. 

The effect of the moon is approximated by adding the mass of the moon 
to the mass of the earth and increasing the initial parking orbit radius 
so that the velocity is the same as in a 100 n. m. earth orbit. Also 
included in the first part are curves showing the variations of AV with 
respect to transfer time and typical two and three-impulse primer vector 
histories. The experience gained in solving the 3-body trajectory 
optimization problems has been most valuable in the subsequent develop- 
ment of the 4-body trajectory optimization program. 

The second part of this report deals with the development of a 
comprehensive program to compute fuel optimal 4-body trajectories 
between the earth and some point in the Sun-Earth-Moon system. The 
moon is treated as a separate entity. The basic building blocks of the 

-program-a-re-the-integ-raton and-the-iterator_s. The manner i n which 

these building blocks are connected depends on the selection of the 
dependent and independent variables. The integrator uses the Stumpff- 
Weiss method^* ^ to extrapolate the state vector and to compute the 
state transition matrix. An important feature is that the cost and 
constraint gradients can be computed analytically in terms of the terminal 
state and the state transition matrix. This method does not require the 
switching of coordinates and generates its own ephemerides. The 
iterators solve the boundary value problems to satisfy terminal conditions 
or to optimize AV with or without terminal constraints. 

The generation of a 4-body fuel optimal trajectory is considerably 
more difficult and time-consuming than in the 2 -body system mainly 
because of the increased difficulty in solving the Lambert problem. In 
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a 4-body 2 -impulse transfer it is generally necessary to go through a 
search process to obtain an initial estimate of the required velocity. If 
the initial estimate is reasonably good so that the terminal miss is small, 
an iterative solution of the boundary value problem will converge to the 
required velocity. If the boundary conditions and/or the transfer time 
are changed, the required velocity for the perturbed trajectory may be 
obtained by iterating on the Lambert solution for the reference trajectory 
provided that the changes are kept small. In general, the solution of the 
Lambert problem will require several iterations, each involving a costly 
function evaluation (the extrapolation of the state vector and the state 
transition matrix). On the other hand, the computation of the gradient in 
terms of the terminal state and the state transition matrix is trivial. 

A number of iterators which have been successful in solving 2 -body 
trajectory optimization problems perform poorly in the 4-body system. 
They either require too many function evaluations or just fail to converge. 
The selection of iterators is thus of major importance. 

Some reduction in computer time is possible by a judicious choice 
of the independent variables. For instance, the classical choice of the 
independent variables to optimize AV in a 2 -body 3 -impulse transfer 
is the position and time of the interior impulse' ' ' . The gradient of 

AV with respect to the independent variables may be expressed in terms 
of the time derivative of the primer vector. Since the solution of the 
2 -body Lambert problem is a single step process, there is no inner loop 
of importance. When this approach is applied to a 4-body problem, it 
would require the solution of two difficult inner loop Lambert problems 
to satisfy constraints at two places and an outer loop to optimize AV. 

The 4-body problem is highly non-linear in that the inner loops will fail 
to converge unless the changes in the interior impulse position and time 
as generated by the outer loop are heavily constrained. As a result, 
the progress to a converged solution tends to be very slow. 

In view of the fact that 'the reference and the perturbed, trajectories 
are required to satisfy the same boundary conditions, there are only 
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four degrees of freedom. Thus, a better approach is to iterate on the 
initial required velocity and the time of the interior impulse in the outer 
loop. The effect of this change is that one of the two inner loop Lambert 
problems is eliminated. The gradient of cost with respect to the new 
independent variables may be computed without computing the primer 
vector. This new approach results in a significant saving in computer 
time. To insure convergence the required change in the interior impulse 
position with respect to the reference trajectory for the remaining 
Lambert problem is introduced in increments rather than in one single 
step. After the problem has converged to a solution, the primer vector 
history is computed to determine whether the trajectory is optimal or an 
additional impulse is required. 

The 4-body trajectory optimization program provides the capability to 
compute 2-impulse transfers between the earth and a point in space with 
or without optimization and 3 -impulse fuel-optimal transfers between 
two points in space. The 2-impulse transfer may originate from the 
earth or a point in space. The terminal condition of a point to earth 
transfer may be either of a point to point type (PTP) or of a point to 
radius type (PTR). The initial condition of a transfer from the earth to a 
point in space is of a PTP type in which the initial position is also varied 
to optimize AY. 


This report is concluded, by showing an example of a 4-body 2- 
impulse transfer from the L^ libration point to the earth. The transfer 
time is chosen to be the same as the typical loop type 3 -body trajectory. 
The difference between the 4-body and 3-body trajectories is readily . 
discernable. 
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PART I 

EXAMPLES OF THREE- BODY TRAJECTORIES 
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1. 1 Generation of Three-Body Trajectories 

Both 2-impulse and 3-impulse trajectories in the 3-body system 
have been generated between the libration point on the Sun-Earth line 
and the earth. The motion of the earth is assumed to be circular around 
the sun. The trajectories of the earth and the vehicle are confined to the 
ecliptic. The 2 -impulse transfer originates from the L^ point and 
terminates at the earth with a given radius and zero radial velocity. The 
terminal radius corresponds to a parking orbit which would give the 
same circular velocity as in an 100 n. m. orbit after adding the mass of 
the moon to the mass of the earth. 

The procedure used to generate a typical family of 2 -impulse 
trajectories are as follows. 

1. A magnitude of the required velocity at L^ is selected and the 
direction relative to the Sun-Earth barycentric line is varied. 

The state vector is extrapolated to a prescribed terminal time 
using the Wilson's version of the multi-conic method^ . The 
state transition matrix is computed along with the trajectory. 

A field of trajectories is generated in this manner. An initial 
estimate of the required velocity is obtained from the trajectory 
which has the smallest miss distance. , 

2. T he const raint g rad ient is com pute d in term s of the termina l 

state and the state transition matrix. The Lambert problem 
is then solved by iterating on the initial velocity using the 
Newton- Raphson method until the desired terminal conditions 
are satisfied. 

3. The terminal time is varied and the Lambert problem is 
re-solved using the required velocity of the reference 
trajectory as the initial estimate. A parabola fit is used 
to extrapolate the required velocity from three adjacent 
Lambert solutions. 

4. The variations of AV with respect to flight time is plotted. 

The primer vector history of the two-impulse transfer which 
requires the minimum AV with respect to the flight time is 
examined. 
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5. If the magnitude of the primer vector history exceeds one, a 
three-impulse trajectory is generated using a method similar 
to the one described under Part II. 

I. 2 Examples of Typical 3-Body Trajectories from to Earth 

The 3 -body trajectories may be grouped into four distinct families. 
Typical trajectories of these families are summarized in Table I. 

The first three families of trajectories have corresponding 
trajectories found by previous investigators for moon to Lg transfers 
in the Earth-Moon system. The last family is new. 
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TABLE I 


Typical 3 -Body Trajectories 


Family 

Typical 2-Impulse Trajectories 

AV 

vs 

Time 

Primer 
Vector 
History for 
3 -Impulse 
Trajectory 

. 

Transfer 

Time 

(Day) 

Traj. 

AV 

(m/s) 

Primer 

Vector 

History 

Fast (Lj to 

leading edge 
of earth) 

33. 563 

Fig. 1 

454. 446 

Fig. 2 

Fig. 5 

Note 2 

Fast (L-^ to 

trailing edge 
of earth) 

35. 563 

Fig. 3 

341.350 

Fig. 4 

Fig. 5 

Note 2 

Loop type 

116. 682 

Fig. 6 

272.137 

Fig. 7 

Fig. 8 

Fig. 9 

Double pass 
type 

174. 32 

Fig. 10 

203. 34 

Fig. 11 

Fig. 12 
Note 4 

Note 3 


NOTE: 

1. Mass ratios used are shown in Appendix A. 

2. The 2 -impulse transfer is fuel- optimal. A third impulse is not 

_ _ _ required. 

3. A converged solution to a 3 -impulse transfer has not been obtained. 

4. Extrapolation has not been carried out far enough to establish 
the flight time for minimum AV. 
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II. 1 General 


The basic building blocks of the 4-body trajectory optimization 
program are the integrator and the iterators. The integrator is the 
routine to extrapolate state vector and. compute the state transition 
matrix. The iterators are the routines to solve boundary value problems 
with or without optimization. The types of problems which will be 
encountered may be classified into three groups. 

1. No optimization - Solution of Lambert problem to satisfy terminal 
constraint only. 

2. Optimization only - Example - outer loop of a 3 -impulse transfer 

3. Optimization with constraints - Example - two-impulse transfer 
from the earth to a point in space 

The structure of the program is determined by the choice of the 
dependent and independent variables. 

II. 2 C oordinate System 

All coordinate systems are three dimensional, non-rotating, 
parallel Cartesian systems with their origins located at the center of 
mass of the body. The x-y plane of the sun-centered system lies in the 
ecliptic. The orientation of the axes in the ecliptic is not specified and 
will depend on the source which provides the ephemerides for the initial 
conditions. A possible choice, for example, is to require the x-axis to 
intersect the Vernal Equinox at the beginning of the nearest Besselian 
year. 

II. 3 Integrator 

There are only two methods which are suitable for use in a 4-body 

trajectory optimization program. They are both multi-conic methods. 

( 8 ) 

The first version is due to Wilson v ' and the second version is due to 
Stumpff and Weiss^^. The major difference between the two methods 
is that during each time step the former computes the two-body conics 
in sequence while the latter computes them in a parallel manner. In a 
4-body space the Stumpff-Weiss method has the following advantages: 
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1 . 


At each time step, no logic is required to determine the proper 
sequence of computation of the conics. While this is not a problem 
in a 3-body space, it introduces complications in a 4-body space. 

2. The method generates its own ephemerides for the 4-body space. 

If the ephemerides are read, say, from a JPL tape, it is not 

possible to exclude the presence of the other bodies. 

Input ephemerides are needed as initial conditions only. The 
required inputs are the initial positions and velocities of the earth, moon 
and vehicle with respect to the sun at the initial time and the mass ratios 
of the massive bodies. 

II. 3.1 Stumpff- Weiss Method^^ 

The 4-body space is shown in Figure 13. The position and velocity 
vectors of a body are indicated by the second subscript with respect to 
another body indicated by the first subscript. 

At each time step it is required to compute six two-body conics and 
three approximate perturbation vectors. From these conics and perturba 
tion vectors are obtained six reference trajectories. The state transition 
matrix is computed from the two-body state transition matrices. The 
reference trajectories are in error equal to the remainders which will 
be computed every four time steps to indicate the errors but not used to 
correct the reference trajectories. This is because the computed, state 
transition matrix is valid only for the reference trajectories. If the 
errors are too big, the user may decrease the step size or increase the 
number of steps accordingly. 

At the first time step, the reference trajectories equal the true 
trajectories. At subsequent time step, the reference trajectories are 
used as inputs to generate the conics and the approximate perturbations. 

(13) 

The six conics will be generated using Goodyears routines . 

They are denoted by: 
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lp i, j ] = ([F i, i l [v i, j" 
i = S, E, M 
j = E, M, V 

i t i 

The six reference trajectories are denoted by: 



i 


3 = 
i t 


S, E, M 
E, M, V 
3 


The three approximate perturbations are given by equations below. 


Me 


SV = pig + pi E ^ P SE^ " J P SE (0) ^ + ^ P EV^ " J P EV 


( 0 )} 


pig +^i M ^ P SM ] ' J P SM (0) ^ + ^MV 1 “ J P MV (0) ^ 


( 1 ) 




P SE pig + pi M ^ P SM^ “ J P SM (0) ^ ’ pi M +^ E ^ P EM^ “ J P EM (0) ^ 


P SM " Mg + M e ^ P SE ] “ J P SE (0) } + pi E + pi M ^ [p EM ] * J P EM (0) ^ 


where 

pi^ = mass ratio, i = S, E, M 
h = step size 


p lj( °) 


reference trajectory at beginning of step 


i = S, E, M 
j = E, M, V 
i i 3 
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J 3 hI 3 


°3 l 3 


The reference trajectory is computed from Equation (2). 

Psv " [Psvl + Psv 

PSE = ^SE^ + P SE 
PSM = ^SM 1 + P SM 

PeV “ PsE + Psv 


( 2 ) 


P MV " P SM + P 


SV 


P EM " P SE + P SM 


The state transition matrix is computed from state transition matrices 
of two-body conics in Equation (3)^\ 


cp^(t+h, t) = <p sv (t+h, t) + cp EV (t+h, t) + <p MV (t+h, t) - 2 J 
cp(t+h, t Q ) = <P h (t+h, t) <p(t, t Q ) 


(3) 


The errors or remainders in the reference trajectories propagate according 
to Equation (4). 

* « ' 

R SV = S SV + ^E (s SE + S EV ) + (s SM + S MV ) 


P SE + + ^ s sm “ S TTM^ 


S P E 7 °SE SM EM 


( 4 ) 


R SM (f *S + P M* S SM + ^E (s SE + S EM* 


25 



where 




i = S, E, M 
j = E, M, V 

i t j 


( 5 ) 


Since we have used r . instead of the true r„ in s^, the error in 

_ _* 4 ii. 4 

(s.. - s ..) is h . Then the errors in R is also h and the errors in R 

g 

is h if the remainders are used to correct the reference trajectories. 


The values of R for four time steps are saved and the remainders at 

every fourth step are computed using Stirling's five-point formulation 

( 12 ) 

shown in Equation 6 . 


R (t Q + 4h) 


R (t Q + 4h) 


[64 R (t Q + h) + 24 R (t Q + 2h) 

+ 64 (t Q + 3h) + 14 t (t + 4h) 

h 2 ». 

[T92 R (t + h) + 48 R (t + 2h) 

40 0 O 


( 6 ) 


+ 64 R (t + 3h) 


where t = time at beginning of every 4 time steps. 


II. 4 Structure 

The structure of the trajectory optimization program is strongly 
influenced by the choice of the dependent and independent variables. A 
particular choice determines how many integrators and iterators are 
required and how they are connected together. 

In a two-impulse transfer, the logical choice is the required 
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velocity at the initial time. If the transfer time is fixed and the transfer 
originates from the earth, the initial position vector must also be 
iterated. The terminal constraint may be of a point to point (PTP) type 
or a point to radius (PTR) type. The PTR terminal consists of three 
components: the desired radial distance at earth, the desired radial 
velocity (usually zero), and the desired orbital inclination angle with 
respect to the ecliptic. The PTP terminal is the desired position vector. 
The two-impulse transfer is shown in Fig. 18 and 19. 

In a three-impulse transfer, the problem is to perturb a reference 
trajectory in such a way that cost is reduced.- Both the reference 
trajectory and the perturbed trajectory must satisfy the same PTP 
terminals. 

(2 3) 

In the classical method ’ , the independent variables are the 

position and time of the interior impulse t^). ft requires two 

inner loops to solve two Lambert problems to satisfy terminal conditions 

at t and t P and an outer loop to change (F , t ) to reduce AV. This 
m i mm 

method was developed for two-body problems where the Lambert problem 

can be solved in one step. In the four-body problem the solution of the 

Lambert problem must be done in multiple steps and is very time 

consuming. In adapting the Optimum Multi-Impulse Rendezvous Program 
(4) 

(OMIR) to a four-body problem it was found that most of the computer 
time was spent in extrapolating the state vector and solving the Lambert 
problem in the one -dimensional search. Furthermore, the four-body 
transfer trajectory is very sensitive to changes in F m and t . The 
Lambert problems will converge only when AF m and At m are kept very 
small. The net effect is that the reduction in AV is small from one 
iteration to another. 

The family of the three-impulse transfer trajectories has four 
degrees of freedom. It is then possible to change the outer loop indepen- 
dent variables from (r , t ) to (v , t ). In this choice, F is deter- 

mm o m + m 

mined by v q and t m . It is only necessary to iterate on to satisfy 
the terminal constraints. In other words, there is only one Lambert 
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Stop 


Figure 14 TWO-IMPULSE TRANSFER 
WITHOUT OPTIMIZATION 




Stop 


Figure 15 TWO-IMPULSE TRANSFER 
WITH OPTIMIZATION 
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t . V . V 

m’ o’ m 



Newton- 

Eaphson 

Step 


r f , v f 


Integrator 

Advance 
State from 
t to t- 


Constraint 

Satisfied 


Stop 

Criterion 

Satisfied 


Compute 

Primer 

Vector 

History 


Inner Loop 
Terminal 
Constraint 
PTP only 


Stop 


Figure 16 THREE-IMPULSE TRANSFER 
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problem to be solved. Tests have shown that it is possible to reduce 
computer time by 50%. 

The proposed program structures for 2- and 3 -impulse transfers 
are shown in Figs. 14, 15 and 16. 

II. 4.1 Two-Impulse Transfer from Earth to a Point in Space 

It is assumed that the initial position of the vehicle lies on a 
circular earth parking orbit of known radius r Q and is specified by two 
angles, a and p. (Fig. 17) The initial required velocity is assumed to 
be normal to the position vector and makes an angle y with an intermediate 
plane. The X^Y^Z^, coordinates are inertial non-rotating coordinates. 
The three angles and the magnitude of the required velocity, v , are to 
be chosen so that the cost of transfer from the earth to a specified point 
(P) in a specified time is optimized (Fig. 18). 


II. 4. 1.1 Estimate of Initial Required Velocity at Earth 

An estimate of the initial required velocity at earth for a 2-impulse 
transfer to a point in space may be generated as follows. Let 


EV 


position vector in X^Y^Z^, coordinates 


r Q = magnitude of position vector 
Rotation matrices (from X 'Y 'Z ' to X^Y^Z^: 



/ cos a 

-sin 

a 0^ 


sin a 

cos 

a. o 


1° 

0 

1 


\ 




1 cos P 

0 

-sin fi \ 

% ° 

0 

1 

0 


i sin /3 

0 

COS P 


\ / 
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Then 




0 

0 

= 

[ 

0 

I 

cos y 

-sin y 


1° 

sin y 

cos y 



/r \ 



EV ~ E a R j3 1 0 

i 

0 


cos a cos £ 

sin a cos j3 

sin j3 


EV 


\ r 

\ o 

R R fl R„ v 
a j3 y I o 

0 I 


v q (-sin a cos y - cos a sin j3 sin y ) ^ 
v q (cos a cos y - sin a sin £ sin y) 
v q cos j3 sin y j 


r o ~ r SV = r SE + r EV 


V = V = V + V 

o SV SE EV 


where r Q and v q are position and velocity vectors relative to the sun in 
X-Y-Z coordinates. 


For a specified transfer time and magnitude of the initial impulse 
the initial required velocity is expressed in terms of the 3 angles. By 
varying a, j3 and y in small increments a field of trajectories is generated. 
A crude search is made first for the required velocity which gives the 
minimum miss distance. This is followed by a refined search around 
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the above minimum. 


II. 4. 1. 2 Cost and Constraint Gradients 
Independent variables: 


x = ( v qJ a, 0, y) 


State transition matrix: 


t^) 


*11 

*12 \ 

*21 

^22 i 


Terminal condition: 

7 i - F f d 

Variational equations: 

“f ' *11 SF o + *Pl2 67 o 


6v f = *21 Sr o + *22 6v o 


The variations of the initial position and velocity are related to the 
changes in the independent variables by 


6r, = 


6r 

dr 

= — dx 

o 

dx 

6v 

dv 

= — dx 

o 

dx 

dr 

dv 

o 

+ 

dx 


■ dm 


dx 


dx 
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6v f = |_cp 


3r 


3v 


21 — + ^22 _ 
ax ax 


u V -1 

o 


av 


dx = 


f dx 


3x 


aF 


ax 


0 

0 

0 


-r sin a cos 0 
r Q cos a cos 0 


-r . cos a sin 8 
o r 

-r sin a sin 8 
o r 

r cos 8 
o 


O' 

0 

0 




dv 


3x 


*-sin a cos y - cos a sin 0 sin y 
cos a cos y - sin a sin j3 sin y 
cos /3 sin y 


v q (-cos a cos y + sin a sin £ sin y) 
vo (-sin a cos' y- cos a sin 0 sin y) 
0 


■v cos a cos j9 sin y 
o 

•v sin a cos j3 sin y 

•v sin 8 sin v 
o r ' 


v q (sin a sin y - cos a sin 0 cos y)* 

v q (-cos a sin y - sin a sin /3 cos y) 

v cos B cos y 
o 


Cost gradient: 



Av = 

l y fd " 

v fl 



= 

kv f | 



& |Av f | 

a- T 

. “ Av f 

= 

— T 
-Av f 

r ar av 

f (Q __2 + to 

dx 

|Av f | 

Sx 

|Av f | 

L^21 -- ^22 ~ 

ax 3 x 

3 Av 

a|£v f | 





3x 


3x 


Note that only the Av at the terminal point is considered. 


35 



Constraint gradient: 


' b , Sr f 3r o . rt 3v o 

-SL - — (r f - r fd ) = ■ — = + <Pi2 ~ 

&x ax 1 ia ax 11 ax ax 


11.4,2 Two-Impulse Transfer from a Point in Space to Earth 


The 2 -impulse transfer from a point in space to the earth is shown 
in Fig. 19. An estimate of the initial velocity may be generated in a 
manner similar to the transfer from the earth to a point in space. 

Independent variables.- 

. x = v Q 
State Transition matrix: 


cp(t f , t Q ) = 


^11 ^12 


^21 ^22 


Variational equations: 


5r f = ^12 


<"t = * 22 6v o 


It follows. 
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Cost gradient: 


Av = |AvJ + |Av f | = |v q - v Q "| + |v fd - v f 


where 


- )M e ’ _ _ 

v fd ' V— unit «r f x v ( ) x r f ) 


+ V 


Earth 


SE 


v c _, = velocity to earth 
oil 


r fE = rac ^ a l distance of vehicle at earth 

— T — T — 
dAv Av q Av f 5v f 

3V o ' AV ol l A ~fl aV 0 

_ T T 

AVq Ay f 


|AvJ |Av. 


ITT ^22 


Point to Radius Terminal Condition (PTR) 
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fd 


fd 


- r r 


- r. 


cos 1 d ~ cos 1 


Earth 


where 


r f = l r f I 


. - T - , 

r f r f v f' r f 


cos i = h/h 


h = r f x v f 


u z = (0, 0, 1) 
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0 
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y v 

y f 


0 -x f 


then, 
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-5L = r ' x -1 - V ' x _i 

Sv dv dv 

O O O 


= R 


^22 " V ^12 


Sh 

ST 


— 1 T — 
n Sh 

1 T ~ 

dv 


h T 


= IT (R ^22 - V *12> 


u T 


5 cos i _ z 
3 v_ h' 


(h 2 I - hh T ) (R <p 22 - V <p 12 ) 


Point to Point Terminal Condition (PTP) 


0 (t f ) = r f - r.j 


fd 


it = !!l s 


”12 


II. 4. 3 Three-Impulse Transfer 


(2,3) 


A 3 -impulse transfer between 2 points in space is shown in Fig. 20. 
Given a reference trajectory r, the problem is to generate a neighboring 
trajectory T ' satisfying the same boundary conditions at reduced AV. 


II. 4 . 3.1 Classical Three-Impulse Transfer 

In the classical approach, there is an outer loop which iterates on 

7 and t to reduce AV and two inner loops which iterate on "v and 
m m o 

V to satisfy boundary conditions at t m and t^ respectively. The cost 

gradient of the outer loop is expressed in terms of the primer vector 

and its time derivatives. 
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3Av 


at 


m 



+ 

V 

m 
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11.4.3.2 New Three-Impulse Transfer 

The independent variables in the outer loop are changed to and 
t . The state vector is extrapolated to t . There is only one inner 
loop which iterates on v r m + to satisfy terminal condition r^. 

Cost and Cost Gradient 


Av = |AvJ + |Av m | + | A' 


= tv - v | + |v + -v | + |v- , - V. 

1 o o' 1 m m, 1 1 fd f 


3Av _ 

av„ 


Av o , A v . iv_) . Av f av t 

|AT o! |AT J ' a7 0 l AV fl a7 0 



Variational Equations for the first leg 


fir = (p .. fir + <p 10 6v = (D 6v 

m r mo, 11 o r mo, 12 o mo, 12 o 


6v - cp fir + <p „„ 6v = <p 00 6v 

m v mo, 21 o r mo, 22- o v mo, 22 o 


fir = 0 
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m 


dv, 


<P 


mo, 12 


av. 


m 


- <r> 


Sv, 


mo, 22 


Sv 


m 


= 0 


at 


m 


Variational Equations for the second leg 


fir. =0 = o. fir + + o, 6v + 

f r fm, 11 m v fm, 12 m 


f ^fm, 21 ^ r m + ^fm, 22 ^ V m 


6v + 
m 


•Pfm, 12 ‘Vm, 11 6r m 


Since 


dr = fir + v dt = fir + + v + dt 
m m m m m m m 


fir + = 6r ” - (v + - v ) dt 
m m m m m 


^ v m - ‘Vm, 12 ‘Vm, 11 ^ r m ^ v m v m ^ 


^fm, 12 ‘Vm, 11 Vo, 12 ^ v o + ‘Vm, 12 ‘Vm, 11 ^ v m v m ^ ^ 
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dv m + 

~^r~ = " ^fm, 12 ‘Vm, 11 Vo, 12 


Sv + 
m 


St 


Vm, 12 Vm, 11 (v m “ v m ) 


m 


6v f = *>tm, 21 l 6r m" ' <v m + ’ V _) ^m 1 


+ 22 ^ ^fm, 12 ^fm, ll^mo, 12 ^ v o +l ^fm, 12 ^fm, 11 ^ v m v m * 


^fm, 21 ' Pfm, 22 ®fm,12 ' f fm, ll 1 ? mo, 12 6v o 


+ [ - «>fm. 21 + ffm. 22 w'' ^fm, U 1 (v m + ' v m' ) d * 
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3V.r 


-1 


— [ Vm, 21 ‘ Vm, 22 Vm, 12 Vm, 11^ ^mo. 


dv 


12 


Sv, 


at 


” ~ [ Vm, 21 “ Vm, 22 Vm, 12 1 Vm, 11 J ( V “ v m ) d:t 


m 


m 


Note: This derivation does not need the computation of the primer 
vector and its derivative. 

11.4.3.3 An Alternate Form of Cost Gradient from Classical Approach 

(2 3) 

In Classical Approach ’ 

on T: J = iAvJ + lAv^l + |Av f | 
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If we use the relation 


-X T ~ 6r " = X T 6v - X T ” 6v 
m m o o m m 

\ T+ 6 T + = -X T 6v f + X T + 6v + 

m m f f m m 


we return to 


6J = X T 6v + \ T (6v + - 6v ~) - X F T 6v f 
o o m m m x x 


which will lead to same gradient as before. However, if we use the 
• *j- ♦ 

expressions for X m and X^ in 


m 


6J - *m + 6F m + - ^m" 6F m' 


we will get a different but equivalent expression for the gradient. This 
is because of X Q and X^ are related through the relation 

X T 6v - X T 6x = C. 


Since 


\n ^mo, 11 \> + %o, 12 X o 


X m ^mo, 21 X o + ^mo, 22 X c 


we have 
X_ 


Vo, 12 1 (X m " ^mo, 11 V 


\n Vo, 21 X o +< ^mo, 22 ^mo, 12 . ^ X m Vo, 11 
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. Thus 
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7~ ^fm, 12 (<p fm, 11 X m " X f )] (v m V m ) dt m 
ot_ 
m 

These cost gradients are different in appearance from the previous 
result. 


The constraint gradient with respect to the required velocity in the 
inner loop of the three-impulse case is the same as shown in the two- 
impulse transfer from a point to earth by iterating on v‘ m + instead of V q . 

II. 4. 3. 4 Solution of Inner Loop Lambert Problem in Increments 

The change in the interior impulse position, dr , between the 
perturbed and reference trajectories due to changes in ~ o and t m is 
usually so large that the "v + on the reference trajectory is a poor 
estimate for the inner loop Lambert problem. This difficulty can be 
alleviated by introducing dr m in increments rather than in one step. 

(Fig. 21). 


to t 


m 


Using v 
to obtain r 


o’ l m 


m 


from the outer loop, we extrapolate the state vector 


dr = r - r ,, 
m mm, old 


Let n = number of increments 


dr 


Ar 


m 


m n 


3 = 0 


Procedure: 


1 . 


r = r . , + Ar 

m, new m, old m 
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2 . 


Av m " <Pfm, 12 ^fm, 11 Ar m 

+ + + 

v = v i t H - Av 

m, new m, old m 

3. Solve inner loop Lambert problem from T to?., using 

+ m, new id ° 

~v as the initial estimate, 

m, new 

4. v ,j = v 

m, old m, new 

r m, old r m, new 

3=3+1 

5. If j < n, go to 1. 

II. 5 Iterators 

The major cost in computer time in a 4-body trajectory optimization 
problem is in function evaluation. Each function evaluation consists of 
the extrapolation of the state vector and the state transition matrix to the 
terminal time in multiple steps. The computation of cost and constraint 
gradients in terms of the terminal state and the state transition matrix is 
trivial. 

The solution of the optimization problem will in general require 
several function evaluations. It is a difficult problem for two main reasons. 

1) . The problem is highly nonlinear. Only small variations from a 

reference trajectory is possible to insure convergence. 

2) . The cost function is often non-unimodal in a search direction. Any 

accelerated gradient method which requires an one -dimensional 
search must use some bracketing procedure to locate even a crude 
minimum in the search direction. A simple bracketing procedure 
will require several function evaluations. As a result, the cost in 
computer time in an one -dimensional search can easily become 
prohibitive. 
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During the early stage of the development considerable effort was 
spent in evaluating iterators which have been used successfully in 
solving 2-body trajectory optimization problems. Some qualitative 
comparisons of their effectiveness in solving the 3-body problems are 
given below. 

II. 5.1 Evaluation of Iterators 

(7) 

1. Jacob son-Oksman Method 

This is an accelerated gradient method, which was successfully 

( 21 ) 

used to optimize three and four-impulse two-body rendezvous trajectories . 
It optimizes a scalar function using cubic fit interpolation in lieu of an 
one-dimensional search. 

If there are constraints to be satisfied, they can be included only 
as a penalty function. The effectiveness of the method is reduced and 
the method cannot correct the general deficiency of the penalty function 
method. 

This method was used to optimize cost in the outer loop of a 
3-body 3-impulse transfer based on the classical formulation. The 
progress was very slow because the interpolated minimum falls very 
close to the starting point. Here, the advantage of not requiring an one- 
dimensional search appeared, to be a disadvantage. 

2. Fletcher- Powell Adaptation of Davidon's First Method^' ^ 

This is a popular accelerated gradient method to optimize a scalar 

function but it requires an one- dimensional search to be effective. The 

algorithm we tested came from a NASA MSC program "Optimum Multi- 

(4) 

Impulse Rendezvous Program (OMIR). It incorporates an one- 
dimensional search method using a Golden Section bracketing procedure 
and cubic fit^^. 

We tested it in the outer loop of a 3-body 3-impulse transfer. It 
was found that the iterator requires an excessive number of function 
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evaluations in the one -dimensional search. About 95% of the computer time 
is spent in this effort and the remaining small portion in updating the 
H-matrix and generating a new direction to proceed. The algorithm 
was satisfactory for the two-body problems where the Lambert problem 
can be solved in a single step and a large number of function evaluations 
is not particularly painful. The computation time can be reduced 
by a change of variables as suggested earlier but the ratio of time spent 
in one -dimensional search and that in generating a new direction will 
remain unchanged. The requirements of a one -dimensional search is 
a definite burden. 

3. Campbell- Moore -Wolf Method (6) 

(17 18} 

This method is equivalent to the Arm strong -Mar qua dt Method^ * ' 

in principle. The algorithm we tested came from the Princeton University 
program "Trajectory Optimization. Program for Comparing Advanced 
Technologies (TOPCATl)"^. It is also a penalty function method to 
handle constraints through internally generated weightings inversely 
proportional to the squares of the allowable tolerances. However, the 
user must supply proper scaling to keep numbers within computable 
range . 

The step size and direction are controlled by an inhibitor, X. When 
X is very small, the method behaves like a Newton- Raphson method to 
satisfy the constraint. When X is very large, the method behaves like a 
gradient method of a weighted scalar product of the constraint violations. 
The value of X is increased when step size is exceeded or if the residual 
grows. It is reduced otherwise. In our test problems, X has a steady 
tendency to grow such that not only the step direction is changed but 
also its magnitude is reduced. Neither characteristic is desirable since 
the user has no further control over the behavior once the mechanism 
for changing X has been programmed and it is difficult to change X in 
an optimal way. 

The method has the inherent difficulties associated with a penalty 
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function method. Furthermore, one would normally desire to take 
gradient steps early in the process and then shift to the Newton- Raphson 
method. This method behaves just in the opposite manner. It is also 
very difficult to provide the proper scaling to be consistent with the 
desired tolerances and computable range. 

4. Modified Fletcher Method 

( 24 ) 

This is another version of Davidon's method proposed by Fletcher 

( 23 ) 

and modified, by Powers' . It still requires a crude one -dimensional 
search. In testing the method in the outer loop of a 3-impulse transfer 
it was found that a crude one -dimensional search as suggested by Powers 
did not work. The reason is that the function being optimized is often 
non-unimodal in the search direction. In order to achieve a cost reduction 
from the starting point it is necessary to use some kind of bracketing 
procedure such as one used in II. 5.1.2. As a result, the number of function 
evaluation required becomes much larger than desired. 

II. 5. 2 Recommended Iterators 

A thorough evaluation of known iterators leads to the following 
selection of iterators. 

1. No optimization 

Newton- Raphson Method with a simple automatic step size control. 

This method will be used in 2 places. 

1. Inner loop of a 3 -impulse transfer 

2. 2 -impulse transfer from a point in space to earth. 

2. Optimization only 

(16 ) 

Davidon’s Second Method' ' This method does not require an one- 

(15) 

dimensional search as in the Davidon's First Method . It has the 
advantage of an accelerated gradient method without the penalty of 
excessive number of function evaluations. It appears to be the most 
efficient method among the ones tested. It will be used in: 

1. The outer loop of a 3 -impulse transfer 

2. A 2 -impulse transfer from the earth to a point in conjunction 
with the gradient projection method. 
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3. Optimization with Constraints 

( 22 ) 

1. Gradient Projection Method This method will be used to 
generate fuel-optimal two-impulse trajectories from the earth to 
a point in space for a fixed transfer time. 

2. Accelerated Gradient Projection Method^* This is a 
accelerated version of the gradient projection method. It is a 
combination of the standard gradient projection method and an 
accelerated gradient method. The Davidon's Second method will 
be used for the accelerated gradient portion of the algorithm. It 
will be used as a replacement of the gradient projection method 
if feasible. 


II. 6 An Example of a Four-Body Two-Impulse Transfer 

An example of a 4-bod.y 2 -impulse trajectory from the libration 
point to the earth has been generated. The transfer time is 116. 682 days 
which is the same as the typical loop type 3 -body trajectory mentioned 
in Part I. This a a 3-dimensional problem. Both the earth and the moon 
have initial position and velocity components normal to the ecliptic. The 
projection of the trajectory with respect to the earth in a plane parallel 
to the ecliptic is shown in Fig. 22. The 3 -body trajectory of the same 
transfer time is shown in Fig. 24. It is evident that the 4-body trajectory 
is noticeably different from the 3 -body trajectory. The presence of the 
moon as a separate entity appears to have a very strong effect on the 
trajectory. The distances of the vehicle from the earth and the moon 
is shown as a function of time in Fig. 25. There are three intervals 
during which the vehicle is closer to the moon than the earth. The mass 
ratios used are shown in Appendix A. The initial conditions of the earth 
and the moon were obtained from a MAC data file on JPL Development 
Ephemeris No. 69. 

The primer vector history of this 4-body trajectory is shown in 
Fig. 23. It is a fuel-optimal 2 1 - impulse transfer for the initial conditions 
used. A search for the optimal initial conditions has not been done. 
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APPENDIX A 



1 A. U. = 1. 49597893 x 10 8 KM 
y T = 1. 001098 x 10" 2 

n 

M e+ m = 3.0404322 x 10 -6 
jj> s = 9. 9999696 x 10 _1 
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Based on mass ratios from JPL 
Tech. Report 32-1306. 
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APPENDIX B 
Accelerations . 
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